'''
例题：
    已知敌方 100 个目标的经度、纬度如23_1.csv所示。我方有一个基地,经度和纬度为(70,40)。
假设我方飞机的速度为 1000 公里/小时。我方派一架飞机从基地出发,侦察完敌方所有目标,再返回原来
的基地。在敌方每一目标点的侦察时间不计,求该架飞机所花费的时间(假设我方飞机巡航时间可以充分长)。
我方有一个基地,经度和纬度为(70,40)
（旅行商问题）

分析：
先编号，基地为1,敌方目标2-101,基地最后再回来为102,所以问题转化为有距离矩阵D=(dij)102*102,
从1走遍所有点到102的最短距离
因为题目给的是经纬度，化简为实际距离d=Rarccos(cos(x1-x2)cosy1cosy2+siny1siny2)
python的三角函数参数是弧度制！！

模拟退火：
1 解空间
  S={(p1,p2,...,p102)|1,(2...101循环),102}
2 目标函数
  min sum(dpip(i+1))  i:1~101
3 新解的产生
  (1)2变换法，任选u,v 对称交换u v中间部分
  (2)3变换法，任选u,v,w，将u v中间的部分放到w后面
4 代价函数差
  delta(f)=dnew-dold
5 接受准则
  P=1                   delta(f)<0
    exp(-delta(f)/T)    delta(f)>=0
6 降温
    T=aT (a=0.999)
7 结束条件
    终止温度 e=10e-30
'''
import pandas as pd
import numpy as np
from math import *
import random
import time
import matplotlib.pyplot as plt


# 手写了一个类似matlab的randperm函数产生随机的全排列
# 用来模拟蒙特卡罗产生搜索初始值
def randPerm(a, b):
    vis = [False for i in range(a, b + 1)]
    res = [0 for i in range(a, b + 1)]
    random.seed(time.time())
    for j in range(0, b - a + 1):
        while (res[j] == 0):
            temp = random.randint(a, b)
            if (vis[temp - a] == False):
                res[j] = temp
                vis[temp - a] = True
            else:
                continue
    return res


# 然而后来发现python还是有这个函数的
list = random.sample(range(2, 102), 100)
# print(list)

def swapUandV(c1, c2, S0):
    St = S0[0:c1]
    St.extend(S0[c2:c1 - 1:-1])
    St.extend(S0[c2 + 1:102])
    return St


data = pd.read_csv('12_1_2.csv')
R = 6371
D = np.zeros((102, 102))
# 填充D矩阵
for i in range(1, 103):
    for j in range(1, 103):
        if i == j:
            continue
        x1 = data['longitude'][i - 1]
        y1 = data['latitude'][i - 1]
        x2 = data['longitude'][j - 1]
        y2 = data['latitude'][j - 1]
        D[i-1][j-1] = R * acos(cos(radians(x1 - x2)) * cos(radians(y1)) * cos(radians(y2)) + sin(radians(y1)) * sin(radians(y2)))

# 蒙特卡罗产生随机初始值
S0 = []
Sum = inf
for i in range(0, 1000):
    S = [1]
    S.extend(random.sample(range(2, 102), 100))
    S.append(102)
    temp = 0
    for j in range(0, 101):
        temp = temp + D[S[j]-1][S[j + 1]-1]
    if temp < Sum:
        S0 = S
        # 这里产生一个Sum，即以当前初始值S0为路线，一直加起来的最大和
        Sum = temp
# 模拟退火
e = 1e-30
at = 0.999
T = 1
# 最大迭代次数
L = 500000
for i in range(0, L):
    # 随机产生u v 2~101
    c = random.sample(range(1, 101), 2)
    c.sort()
    c1 = c[0]
    c2 = c[1]
    # 计算代价函数值
    df = D[S0[c1 - 1]-1][S0[c2]-1] + D[S0[c1]-1][S0[c2 + 1]-1] - D[S0[c1 - 1]-1][S0[c1]-1] - D[S0[c2]-1][S0[c2 + 1]-1]
    # 接受准则
    if df < 0:
        S0 = swapUandV(c1, c2, S0)
        # 相当于减去这个差值，因为Sum是前面算出来的初始值的对应函数值
        Sum = Sum + df
    else:
        if (exp(-df / T) > random.random()):
            S0 = swapUandV(c1, c2, S0)
            Sum = Sum + df
    T = at * T
    if (T < e):
        break
print(S0)
print(Sum)

# 画图
x=[]
y=[]
for i in S0:
   x.append(data['longitude'][i-1])
   y.append(data['latitude'][i - 1])

plt.xlim(0,70)
plt.ylim(0,40)
plt.plot(x,y,color='r', linewidth=1, alpha=0.6,marker='*')
plt.show()

# 搞了一天终于搞定了，基本接近书上的答案44了 恶补了很多东西